Hakalau Isoscape

if (!require('knitr')) install.packages('knitr'); library('knitr')
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')

# Load in packages
if (!require("pacman")) install.packages("pacman"); library(pacman) # for rapid install if not in library

# load packages
pacman::p_load("RgoogleMaps", "SDMTools", "rgdal", "ggmap", "gridExtra", "automap", "cowplot","sp", "gstat", "plotrix", "dplyr", "ggplot2", "scales","magrittr","plyr","effects", "ggpubr", "sjstats", "RColorBrewer")

Background

Hakalau National Forest Wildlife Refuge Isoscapes
Data collected 20-21 August 2019 in Hakalau Forest, Hawai’i Island. Soil and plant samples collected from afforested/restored plots after a century of deforestation from cattle grazing (AK sites) and remnant mixed koa-ohia forests (RK sites). Target engineers of forest habitats are Acacia koa (Hawaiian name: Koa), an indigenous hardwood species that forms native forest canopies and associates with nitrogen fixing bacteria. Understory species include Rubus argutus – an invasive species (common name: sawtooth blackberry) – and an endemic species Rubus hawaiiensis (common name: ’Ākala); both are members of the Roseaceae family. Noteably, R. argutus was found more in the remnant plot (i.e., RK) and the indigenous R. hawaiiensis was common in the restored plot (i.e, AK).

Soil samples were collected in a spatially explicit fashion (every 4 or 5m) with a 20 x 35m superplot, consisting of thirty-five 4 x 5m subplots. This data was used to create a δ¹³C and δ¹⁵N stable isotope landscape, hereafter ‘isoscape’. Koa and Rubus spp. were collected within each 4 x 5m sublot.

Plot densities were…
AK: Koa (18 trees) Rubus (14 plants) per 20 x 35m plot= 0.026 Koa m^-2, 0.020 Rubus m^-2 RK: Koa (10 trees) Rubus (28 plants) per 20 x 35m plot= 0.014 Koa m^-2, 0.040 Rubus m^-2

Site Maps

# load data
Hak.gps<-read.csv("data/Hakalau.gps.csv")
API<-read.csv("data/API_key.csv")
API.key<-API[1,1]

ggplot(Hak.gps, aes(x = longitude, y = latitude)) +
  coord_quickmap() +
  geom_point()

Figure 1. Site map of (a) the Hawaiian Island archipelago, (b) sampling plot in Hakalau Forest Refuge, and (c) plot layouts, highlighting Acacia koa density.

######## using ggmap
register_google(key=API.key)

########
#Hawaiian islands Map
########
hi_map<-get_map(location=c(-157,20.5),zoom=7,maptype="satellite",color="color")
hi_map_for_man <- ggmap(hi_map) +
  geom_point(aes(x = -155.320, y = 19.83), pch=23,colour="black",fill="red", size = 2, stroke=0.5) +
  xlab("") + ylab("Latitude") +
  scale_y_continuous(limits=c(18.8, 22.2))+
  scale_x_continuous(limits=c(-160.2, -154)) +
  theme(axis.text=element_text(colour="black",size=5),
        axis.title=element_text(colour="black",size=8),
        plot.margin = unit(c(0, 0.5, 0, 0.7), "cm")) +
  ggsn::scalebar(x.min=-160.2, x.max=-154, y.min=18.8, y.max=22.2, dist=50, dist_unit="km", transform=TRUE,
                 st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84",st.color="white", border.size=0.5)


##########
# site map plot showing distance from plots
##########
Hak.rev=c(-155.317, 19.82158002)

map2<-get_map(Hak.rev, 
                      zoom=16, 
                      scale = 2, 
                      maptype= "satellite",
                      source="google", extent= "device", legend="topright")

## site single point only
site.only<-Hak.gps[c(1,5),]
Hakalau.site<-
  ggmap(map2)+
  geom_point(aes(x=longitude, y=latitude), data=site.only, alpha=0.5, color="white", size=4)+
  labs(x="", y="") +
  scale_y_continuous(limits=c(19.81785, 19.8242))+
  scale_x_continuous(limits=c(-155.3223, -155.311)) +
  theme(text = element_text(size=6),
       plot.margin = unit(c(0.2, 0.5, 0.2, 0.2), "cm")) +
  annotate("text", x=-155.3187, y=19.8225, label= "AK", size=3, col="white") +
  annotate("text", x=-155.31478, y=19.8224, label= "RK", size=3, col="white") +
ggsn::scalebar(x.min=-155.314, x.max=-155.311, y.min=19.818, y.max=19.824, dist=100, 
               dist_unit="m", transform=TRUE,
              st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84",st.color="white", border.size=0.5)
Hakalau.site

Figure 1. Site map of (a) the Hawaiian Island archipelago, (b) sampling plot in Hakalau Forest Refuge, and (c) plot layouts, highlighting Acacia koa density.

##########
# Site plots: AK or RK only
##########
AK<- Hak.gps[(Hak.gps$Site=="AK"),]
RK<- Hak.gps[(Hak.gps$Site=="RK"),]

####### AK map
AK.gps<-c(-155.3189, lat=19.8219)
mapAK<-get_map(AK.gps,
              zoom=19, 
              scale = 2, 
              maptype= "satellite",
              source="google", extent= "device", legend="topright")

AK.map<-ggmap(mapAK, group=type)+
  geom_point(aes(x=longitude, y=latitude, shape=type, color=type), data=AK, alpha=0.5, size=3)+
  scale_y_continuous(limits=c(19.82165, 19.82225))+
  scale_x_continuous(limits=c(-155.3194, -155.3183)) +
  theme(legend.position="top",
    legend.text = element_text(color = "white"),
    legend.title = element_text(color = "white"),
    legend.key = element_rect(fill = "white"),
    plot.margin = unit(c(0, 0.5, 0, 0), "cm"),
    text = element_text(size=8)) +
  guides(colour= guide_legend(override.aes = list(color = "white"))) +
  scale_color_manual(values=c('red','mediumseagreen'))+
  labs(x="Longitude", y="Latitude") +
  annotate("text", x=-155.31925, y=19.82225, label= "Restored Forest (AK)", size=3, col="white")+
  ggsn::scalebar(x.min=-155.3183, x.max=-155.3187, y.min=19.82165, y.max=19.82225, dist=10, 
               dist_unit="m", transform=TRUE,
              st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84", st.color="white", border.size=0.5)


####### RK map
RK.edit<-RK[-12,] # remove "proximate" Koa
RK.edit[9,3]<-"koa" # rename to be "koa"

RK.gps<-c(-155.315, 19.8216)
mapRK<-get_map(RK.gps,
               zoom=19, 
               scale = 2, 
               maptype= "satellite",
               source="google", extent= "device", legend="topright")

RK.map<-ggmap(mapRK, group=type)+
  geom_point(aes(x=longitude, y=latitude, shape=type, color=type), data=RK.edit, alpha=0.5, size=3)+
  scale_color_manual(values=c('red','mediumseagreen', "dodgerblue"))+
  scale_shape_manual(values=c(16, 17, 3)) +
  scale_y_continuous(limits=c(19.82120, 19.82195))+
  scale_x_continuous(limits=c(-155.31575, -155.3144)) +
  theme(legend.position="top",
        legend.title=element_blank(), 
        legend.key = element_rect(fill = "white"),
        plot.margin = unit(c(0, 0.5, 0, 0), "cm"), 
        text = element_text(size=8)) +
  labs(x="Longitude", y="") +
  annotate("text", x=-155.3156, y=19.82193, label= "Remnant Forest (RK)", size=3, col="white")+
  ggsn::scalebar(x.min=-155.3148, x.max=-155.3144, y.min=19.8212, y.max=19.8219, dist=10, 
               dist_unit="m", transform=TRUE,
              st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84", st.color="white", border.size=0.5)


site.plots<-plot_grid(hi_map_for_man, Hakalau.site, AK.map, RK.map, 
          labels=c('a', 'b', 'c', 'd'), label_size=8, hjust=-1, vjust= 2, ncol=2, nrow=2)

### export it
pdf(file= "figures/Fig1.sitemaps.pdf", height=7, width=8)
site.plots
dev.off()

## quartz_off_screen 
##                 2

Ecology and Isotopes data

• DBH calculations
  
• data mutating
  
• test scatters of all samples by habitat

Hak.df<-read.csv("data/Hakalau.isotopes.csv")

########## DBH 
DBH.df<- Hak.df %>%
  select(c(DBH..cm.1, DBH..cm.2, DBH..cm.3, DBH..cm.4, DBH..cm.5))

# replace NAs with zeros
DBH.df<- DBH.df %>%
  mutate_each(funs(replace(., which(is.na(.)), 0)))

# if only one trunk, then make new column and report this value
DBH.df <- DBH.df %>%
  mutate(DBH.1trunk = ifelse(DBH..cm.2>"0", 0, DBH..cm.1))

# if multiple trunks, sum of squares and square root all trunks
# and plug in the DBH for single trunks
DBH.df <- DBH.df %>%
  mutate(DBH.Total = ifelse(DBH.1trunk=="0",  
                     sqrt(DBH..cm.1^2 + DBH..cm.2^2 + DBH..cm.3^2 + DBH..cm.4^2 + DBH..cm.5^2), DBH.1trunk))

# replace zeros with NA now that total DBH has been calculated
DBH.df$DBH.Total<-as.numeric(ifelse(DBH.df$DBH.Total=="0", "NA", DBH.df$DBH.Total))

############ Canopy diameter
# use area of an elipse, a x b x π, where a and b are major and minor radius
Can.df<-Hak.df %>%
  select(c(canopy.0..m, canopy.90..m, canopy.180..m, canopy.270..m))

# Major radius is the average of 0 and 180 degrees
Can.df<- Can.df %>% 
  mutate(Maj.rad = rowMeans(cbind(canopy.0..m, canopy.180..m), na.rm=T))

# Minor radius is the average of 90 and 270 degrees
Can.df<- Can.df %>% 
  mutate(Min.rad = rowMeans(cbind(canopy.90..m, canopy.270..m), na.rm=T))

# Elipse area
Can.df<- Can.df %>% 
  mutate(Canopy.area = ifelse(Maj.rad>0, Maj.rad*Min.rad*3.14, "NA"))

########################
# combine DBH.Total and Canopy.area with isotope dataframe and reduce columns
Hak.df$DBH.Total<-DBH.df$DBH.Total
Hak.df$Canopy.area<-Can.df$Canopy.area

# Hak.data is the full dataframe
Hak.data<- Hak.df %>%
  select(c(Plot, Position.point, Sample, DBH.Total, Canopy.area, d15N, d13C, N..percent, Total.N..mmol.gdw, C..percent, Total.C..mmol.gdw, C.N))

########################
########################
########################

# test plots

### Scatter of d13C vs. d15N
iso.scatter<-ggplot(data=Hak.data, aes(x=d13C, y=d15N, color=Plot, shape=Sample))+
  geom_point(cex=1.5) +
  xlab(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) 

### Scatter of total C vs. total N
total.nut<-ggplot(data=Hak.data, aes(x=Total.C..mmol.gdw, y=Total.N..mmol.gdw, color=Plot, shape=Sample))+
  geom_point(cex=1.5) +
  xlab("total Carbon") +
  ylab("total Nitrogen")

#  figures together
plot_grid(iso.scatter, total.nut,
     labels=c('a', 'b'), label_size=10, hjust=-1, vjust= 3, ncol=2, nrow=1)

Rubus and Koa relationship

• plotting the relationship between Rubus and Koa d15N
  
• took average of 3 closest Koa trees to each Rubus plant
  
• plot Rubus d15N against Koa d15N (mean of 3 trees)

Ndfa<-read.csv("data/Ndfa.Rub.csv")

# rename the d15N so we know it is for rubus
Ndfa$d15N.Rub<-Ndfa$d15N
Ndfa$d15N.Koa<-Ndfa$mean.d15N.Ndfa

####################
#################### relationship between koa and rubus d15N?
# Yes, but a bit stronger in the RK site
mod.AK<-lm(Ndfa[(Ndfa$Plot=="AK"),]$d15N.Rub~Ndfa[(Ndfa$Plot=="AK"),]$d15N.Koa); anova(mod.AK)

## Analysis of Variance Table
## 
## Response: Ndfa[(Ndfa$Plot == "AK"), ]$d15N.Rub
##                                      Df Sum Sq Mean Sq F value  Pr(>F)  
## Ndfa[(Ndfa$Plot == "AK"), ]$d15N.Koa  1 1.9786 1.97856  4.4696 0.05612 .
## Residuals                            12 5.3120 0.44267                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

int.rubAK<-coef(summary(mod.AK))[1] #intercept
slop.rubAK<-coef(summary(mod.AK))[2] #slope
                                  
mod.RK<-lm(Ndfa[(Ndfa$Plot=="RK"),]$d15N.Rub~Ndfa[(Ndfa$Plot=="RK"),]$d15N.Koa); anova(mod.RK)

## Analysis of Variance Table
## 
## Response: Ndfa[(Ndfa$Plot == "RK"), ]$d15N.Rub
##                                      Df Sum Sq Mean Sq F value   Pr(>F)   
## Ndfa[(Ndfa$Plot == "RK"), ]$d15N.Koa  1  5.503  5.5030  8.9346 0.006043 **
## Residuals                            26 16.014  0.6159                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

int.RKoa<-coef(summary(mod.RK))[1] #intercept
slop.RKoa<-coef(summary(mod.RK))[2] #slope

#### plot
BW.back<-theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black", size=0.2))

Rub.Koa.d15N.plot<-ggplot(data=Ndfa, aes(x=d15N.Koa, y=d15N.Rub, color=Plot))+
  geom_point(aes(color=Plot), alpha=.8) + 
  scale_color_manual(values=c('gray30', "gray70")) +
  coord_equal() + theme_bw() +
  xlab(expression(paste(italic("Acacia koa "), delta^{15}, N, " (\u2030)"), sep=""))+
  ylab(expression(paste(italic("Rubus"), " spp. ", delta^{15}, N, " (\u2030)"), sep="")) +
  guides(colour=guide_legend(override.aes = list(stroke=1, fill=NA, alpha=1))) +
  geom_smooth(method = "lm", alpha = .15, size=0.6) +
  theme(legend.title = element_blank(),
        axis.text=element_text(size=6),
        axis.title=element_text(size=10))+
  annotate(geom="text", x=2.7, y=0.8, label=expression(paste("AK ", italic("p"), "=0.056")), size=2) +
  annotate(geom="text", x=2.7, y=0.6, label=expression(paste("RK ", italic("p="), bold("0.006"))), size=2) +
  BW.back

Rub.Koa.d15N.plot

Figure 2a. Relationship between d15N-Koa and 15N-Rubus at two sampling sites.

dev.copy(pdf, "figures/Fig 3.Koa.Rub.regr.pdf", width=4, height=4, encod="MacRoman")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

DBH, Canopy d15N models

• plot against d15N values with model curves
• scatters with models
• d15N differ by PLOT (AK vs RK), no other effects or relationship with DBH or Canopy

# ** No relationship between d15N and DBH or Canopy at plot level**
# ** Effects at PLOTs on d15N (lower in AK) **


######## plot of d15N and DBH, and d15N and Canopy area

# all samples
par(mfrow=c(1,2), mar=c(4,5,1,1))

#### DBH
mod.DBH<-lm(d15N~DBH.Total*Plot, data=Hak.data) # no interaction effect

mod.DBH<-lm(d15N~DBH.Total+Plot, data=Hak.data)
int<-coef(summary(mod.DBH))[1]
DBH.slope<-coef(summary(mod.DBH))["DBH.Total",1] # DBH (slope for all figures)
RK.int<-coef(summary(mod.DBH))["PlotRK",1] # RK, AK = 0


# no significant relationship to d15N but a trend for DECREASE with DBH
plot(d15N~DBH.Total, data=Hak.data,
    col=c("coral","dodgerblue")[as.factor(Plot)],
    pch=c(1,2)[as.factor(Plot)],
    ylab=expression(paste(delta^{15}, N, " (\u2030, air)")),
    xlab="DBH (cm)",
    ylim=c(0,5), xlim=c(0,60))
    ablineclip(int, DBH.slope, col="coral", lwd=2, 
               x1 = min(Hak.data$DBH.Total[Hak.data$Plot=="AK"], na.rm=T), 
               x2 = max(Hak.data$DBH.Total[Hak.data$Plot=="AK"], na.rm=T)) # AK model
    ablineclip((int+ RK.int), DBH.slope, col="dodgerblue", lwd=2, 
               x1 = min(Hak.data$DBH.Total[Hak.data$Plot=="RK"], na.rm=T), 
               x2 = max(Hak.data$DBH.Total[Hak.data$Plot=="RK"], na.rm=T)) # RK model

    
#### Canopy
mod.Can<-lm(d15N~Canopy.area*Plot, data=Hak.data) # no interaction

mod.Can<-lm(d15N~Canopy.area+Plot, data=Hak.data)
int<-coef(summary(mod.Can))[1]
Can.slope<-coef(summary(mod.Can))["Canopy.area",1] # DBH (slope for all figures)
RK.int<-coef(summary(mod.Can))["PlotRK",1] # RK, AK = 0

# no significant relationship to d15N but a trend for INCREASE with canopy area
plot(d15N~Canopy.area, data=Hak.data, yaxt="n",
     col=c("coral","dodgerblue")[as.factor(Plot)],
     pch=c(1,2)[as.factor(Plot)],
     xlab=expression(paste("Canopy (m"^2,")")),
     ylab="", 
     ylim=c(0,5), xlim=c(0,80))
     axis(side=2, labels=F)
     ablineclip(int, Can.slope, col="coral", lwd=2, 
               x1 = min(Hak.data$Canopy.area[Hak.data$Plot=="AK"], na.rm=T), 
               x2 = max(Hak.data$Canopy.area[Hak.data$Plot=="AK"], na.rm=T)) # AK model
     ablineclip((int+ RK.int), Can.slope, col="dodgerblue", lwd=2, 
               x1 = min(Hak.data$Canopy.area[Hak.data$Plot=="RK"], na.rm=T), 
               x2 = max(Hak.data$Canopy.area[Hak.data$Plot=="RK"], na.rm=T)) # RK model
     legend("topright", legend=levels(Hak.data$Plot), box.lty=0, bg="transparent", lty=1, pch=c(1,2), 
            col=c("coral", "dodgerblue"), cex=1)

dev.copy(pdf, "figures/d15N.DBH.Canopy.pdf", width=7, height=5, encod="MacRoman")
dev.off()

With Plot being most important effect, divide DBH and Canopy to be plot mean +/- SE
Supplemental Figure 1. with afforested and remnant forests
- DBH and Canopy figure as bar chart and means - using Mann-Whitney nonparametric tests: significant plot effects - greater DBH (AK), trend for greater Canopy (RK)

## Models
mod.DBH<-lm(DBH.Total~Plot, data=Hak.data); anova(mod.DBH) # signif

## Analysis of Variance Table
## 
## Response: DBH.Total
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Plot       1 1290.6 1290.58   8.369 0.00762 **
## Residuals 26 4009.5  154.21                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod.Can<-lm(Canopy.area~Plot, data=Hak.data); anova(mod.Can) # not signif, a bit wobbly in assumptions

## Analysis of Variance Table
## 
## Response: Canopy.area
##           Df Sum Sq Mean Sq F value Pr(>F)
## Plot       1  654.2  654.16  2.0862 0.1606
## Residuals 26 8152.6  313.56

# diagnostic plots
for(i in c(4:5)){
  Y<-Hak.data[,i]
  full<-lm(Y~Plot, data=Hak.data, na.action=na.exclude)
  R <- resid(full) #save glm residuals
  op<-par(mfrow = c(2,3), mar=c(5,4,1,2), pty="sq")
  plot(full, add.smooth = FALSE, which=1)
  QQ <- qqnorm(R, main = colnames(Hak.data)[i]) 
  QQline <- qqline(R)
  hist(R, xlab="Residuals", main = colnames(Hak.data)[i])
  plot(Hak.data$Plot, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
}

# Mann-Whitney U-Test
mwu(Hak.data, DBH.Total, Plot) # signif

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 18) | 2 = RK (n = 10):
##   U = 315.000, W = 144.000, p = 0.010, Z = 2.589
##   effect-size r =  0.489
##    rank-mean(1) = 17.50
##    rank-mean(2) =  9.10

mwu(Hak.data, Canopy.area, Plot) # not

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 18) | 2 = RK (n = 10):
##   U = 252.000, W = 81.000, p = 0.666, Z = -0.432
##   effect-size r =  0.082
##    rank-mean(1) = 14.00
##    rank-mean(2) = 15.40

#### plots
DBH.m<-aggregate(DBH.Total~Plot, data=Hak.data, mean)
Can.m<-aggregate(Canopy.area~Plot, data=Hak.data, mean)

DBH.SE<-aggregate(DBH.Total~Plot, data=Hak.data, std.error)
Can.SE<-aggregate(Canopy.area~Plot, data=Hak.data, std.error)

DBH.n<-aggregate(DBH.Total~Plot, data=Hak.data, length)
Can.n<-aggregate(Canopy.area~Plot, data=Hak.data, length)

mean.eco<-cbind(DBH.m, Can.m[2], DBH.SE[2], Can.SE[2]); colnames(mean.eco)<-c("Plot", "DBH", "Can", "D.SE", "C.SE")

pd <- position_dodge(0.7) #offset for error bars

BW.back2<-theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black", size=0.4))

### DBH plot of mean/SE by plot
DBH.plot<-ggplot(data=mean.eco, aes(x=Plot, y=DBH))+
  geom_bar(stat="identity", position = pd, width=0.7, fill="gray70")+
  geom_errorbar(aes(ymin=DBH-D.SE, ymax=DBH+D.SE), size=.5, width=0, position=pd) +
  xlab("Plots") +
  ylab("dbh (cm)") + 
  theme(legend.title = element_blank(),
        axis.text=element_text(size=6),
        axis.title=element_text(size=8))+
  annotate(geom="text", x=1.5, y=35, label=expression(paste(bold("*"))), size=4) +
  BW.back2

### Canopy plot of mean/SE by plot
Canopy.plot<-ggplot(data=mean.eco, aes(x=Plot, y=Can))+
  geom_bar(stat="identity", position = pd, width=0.7, fill="gray70")+
  geom_errorbar(aes(ymin=Can-C.SE, ymax=Can+C.SE), size=.5, width=0, position=pd) +
  xlab("Plots") +
  ylab(expression(paste("Canopy (m"^2,")"))) + 
  theme(legend.title = element_blank(),
        axis.text=element_text(size=6),
        axis.title=element_text(size=8))+
  BW.back2

#  figures together
plot_grid(DBH.plot, Canopy.plot,
     labels=c('a', 'b'), label_size=8, hjust=-1, vjust= 3, ncol=2, nrow=1)

dev.copy(pdf, "figures/Supp.Fig1.DBHCanopy.pdf", width=4, height=3.5, encod="MacRoman")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

Run statistics on the responses (dd15N, 13C, tissue C and N)

# models

## tests for assumptions
for(i in c(6:12)){
  Y<-Hak.data[,i]
  full<-lm(Y~Plot*Sample, data=Hak.data, na.action=na.exclude)
  R <- resid(full) #save glm residuals
 
  op<-par(mfrow = c(2,3), mar=c(5,4,1,2), pty="sq")
  plot(full, add.smooth = FALSE, which=1)
  QQ <- qqnorm(R, main = colnames(Hak.data)[i]) 
  QQline <- qqline(R)
  hist(R, xlab="Residuals", main = colnames(Hak.data)[i])
  plot(Hak.data$Plot, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
  plot(Hak.data$Sample, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
}

########## Separate dataframes
soil<-Hak.data[(Hak.data$Sample=="Soil"),]
Koa<-Hak.data[(Hak.data$Sample=="Koa"),]
Rubus<-Hak.data[(Hak.data$Sample=="RUBARG" | Hak.data$Sample=="RUBHAW"),]
AK.df<-Hak.data[(Hak.data$Plot=="AK"),]; AK.df.sansoil<-AK.df[!(AK.df$Sample=="Soil"),]
RK.df<-Hak.data[(Hak.data$Plot=="RK"),]; RK.df.sansoil<-RK.df[!(RK.df$Sample=="Soil"),]

d15N model

######## Nitrogen isotopes
d15N.soil<-lm(d15N~Plot, data=soil); anova(d15N.soil) # soil signif.

## Analysis of Variance Table
## 
## Response: d15N
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Plot       1   2.996  2.9963   2.479 0.1187
## Residuals 94 113.612  1.2086

plot(allEffects(d15N.soil))

d15N.Koa<-lm(d15N~Plot, data=Koa); anova(d15N.Koa) # Koa signif.

## Analysis of Variance Table
## 
## Response: d15N
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Plot       1  4.7952  4.7952  7.4173 0.01139 *
## Residuals 26 16.8088  0.6465                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d15N.Koa))

d15N.Rubus<-lm(d15N~Plot, data=Rubus); anova(d15N.Rubus) # Rubus signif.

## Analysis of Variance Table
## 
## Response: d15N
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Plot       1  9.9155  9.9155  13.768 0.0006289 ***
## Residuals 40 28.8076  0.7202                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d15N.Rubus))

d15N.AK<-lm(d15N~Sample, data=AK.df.sansoil); anova(d15N.AK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: d15N
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Sample     1  2.9463 2.94632  6.0621 0.01978 *
## Residuals 30 14.5807 0.48602                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d15N.AK)) 

d15N.RK<-lm(d15N~Sample, data=RK.df.sansoil); anova(d15N.RK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: d15N
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Sample     1  4.4682  4.4682  5.1829 0.02886 *
## Residuals 36 31.0357  0.8621                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d15N.RK))

d13C model

######## ######## 
######## Carbon isotopes
d13C.soil<-lm(d13C~Plot, data=soil); anova(d13C.soil) # soil signif.

## Analysis of Variance Table
## 
## Response: d13C
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Plot       1  24.13 24.1302  22.544 7.342e-06 ***
## Residuals 94 100.62  1.0704                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d13C.soil))

d13C.Koa<-lm(d13C~Plot, data=Koa); anova(d13C.Koa) # Koa signif.

## Analysis of Variance Table
## 
## Response: d13C
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Plot       1 12.825 12.8250  21.018 0.0001006 ***
## Residuals 26 15.865  0.6102                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d13C.Koa))

d13C.Rubus<-lm(d13C~Plot, data=Rubus); anova(d13C.Rubus) # Rubus signif.

## Analysis of Variance Table
## 
## Response: d13C
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Plot       1 21.848 21.8484  46.101 3.684e-08 ***
## Residuals 40 18.957  0.4739                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d13C.Rubus))

d13C.AK<-lm(d13C~Sample, data=AK.df.sansoil); anova(d13C.AK) #only koa vs. rubus comparision

## Analysis of Variance Table
## 
## Response: d13C
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sample     1  0.851 0.85100  1.2075 0.2806
## Residuals 30 21.142 0.70474

plot(allEffects(d13C.AK))

d13C.RK<-lm(d13C~Sample, data=RK.df.sansoil); anova(d13C.RK) #only koa vs. rubus comparision

## Analysis of Variance Table
## 
## Response: d13C
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Sample     1  1.4676 1.46758  3.8622 0.05714 .
## Residuals 36 13.6796 0.37999                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(d13C.RK))

Nitrogen content model

####### Nitrogen content
TN.soil<-lm(Total.N..mmol.gdw~Plot, data=soil); anova(TN.soil) # soil signif.

## Analysis of Variance Table
## 
## Response: Total.N..mmol.gdw
##           Df Sum Sq  Mean Sq F value Pr(>F)
## Plot       1 0.1087 0.108676  1.6099 0.2076
## Residuals 94 6.3456 0.067506

plot(allEffects(TN.soil))

TN.Koa<-lm(Total.N..mmol.gdw~Plot, data=Koa); anova(TN.Koa) # Koa signif.

## Analysis of Variance Table
## 
## Response: Total.N..mmol.gdw
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Plot       1 0.03271 0.032711  0.5942 0.4478
## Residuals 26 1.43136 0.055052

plot(allEffects(TN.Koa))

TN.Rubus<-lm(Total.N..mmol.gdw~Plot, data=Rubus); anova(TN.Rubus) # Rubus signif.

## Analysis of Variance Table
## 
## Response: Total.N..mmol.gdw
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## Plot       1 0.2110 0.21100  3.6973 0.06164 .
## Residuals 40 2.2828 0.05707                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(TN.Rubus))

N.AK<-lm(Total.N..mmol.gdw~Sample, data=AK.df.sansoil); anova(N.AK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: Total.N..mmol.gdw
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Sample     1 1.5930  1.5930  36.455 1.254e-06 ***
## Residuals 30 1.3109  0.0437                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(N.AK)) 

N.RK<-lm(Total.N..mmol.gdw~Sample, data=RK.df.sansoil); anova(N.RK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: Total.N..mmol.gdw
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Sample     1 0.38328 0.38328  5.7415 0.02189 *
## Residuals 36 2.40323 0.06676                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(N.RK))

Carbon content model

######## ######## 
######## Carbon content
TC.soil<-lm(Total.C..mmol.gdw~Plot, data=soil); anova(TC.soil) # soil signif.

## Analysis of Variance Table
## 
## Response: Total.C..mmol.gdw
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Plot       1   81.68  81.678  4.9447 0.02857 *
## Residuals 94 1552.71  16.518                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(TC.soil))

TC.Koa<-lm(Total.C..mmol.gdw~Plot, data=Koa); anova(TC.Koa) # Koa signif.

## Analysis of Variance Table
## 
## Response: Total.C..mmol.gdw
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Plot       1  0.0377 0.03768  0.0962 0.7589
## Residuals 26 10.1849 0.39173

plot(allEffects(TC.Koa))

TC.Rubus<-lm(Total.C..mmol.gdw~Plot, data=Rubus); anova(TC.Rubus) # Rubus signif.

## Analysis of Variance Table
## 
## Response: Total.C..mmol.gdw
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## Plot       1  34.074  34.074  8.6486 0.005419 **
## Residuals 40 157.595   3.940                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(TC.Rubus))

C.AK<-lm(Total.C..mmol.gdw~Sample, data=AK.df.sansoil); anova(C.AK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: Total.C..mmol.gdw
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Sample     1 122.081 122.081  473.95 < 2.2e-16 ***
## Residuals 30   7.728   0.258                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(C.AK)) 

C.RK<-lm(Total.C..mmol.gdw~Sample, data=RK.df.sansoil); anova(C.RK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: Total.C..mmol.gdw
##           Df  Sum Sq Mean Sq F value Pr(>F)  
## Sample     1  32.592  32.592  7.3308 0.0103 *
## Residuals 36 160.052   4.446                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(C.RK))

C:N

####### C:Ncontent
CN.soil<-lm(C.N~Plot, data=soil); anova(CN.soil) # not signif.

## Analysis of Variance Table
## 
## Response: C.N
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Plot       1   3.349  3.3488  1.7346  0.191
## Residuals 94 181.479  1.9306

CN.Koa<-lm(C.N~Plot, data=Koa); anova(CN.Koa) # not signif.

## Analysis of Variance Table
## 
## Response: C.N
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Plot       1   2.289  2.2886  0.3456 0.5617
## Residuals 26 172.171  6.6220

CN.Rubus<-lm(C.N~Plot, data=Rubus); anova(CN.Rubus) # not signif

## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value Pr(>F)
## Plot       1   4.27  4.2660  0.5037  0.482
## Residuals 40 338.78  8.4694

CN.AK<-lm(C.N~Sample, data=AK.df.sansoil); anova(CN.AK) #only koa vs. rubus comparision, signif.

## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## Sample     1  75.46  75.460  12.771 0.001214 **
## Residuals 30 177.26   5.909                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(allEffects(CN.AK)) 

CN.RK<-lm(C.N~Sample, data=RK.df.sansoil); anova(CN.RK) #only koa vs. rubus comparision, not signif

## Analysis of Variance Table
## 
## Response: C.N
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sample     1  24.48  24.482  2.6413 0.1128
## Residuals 36 333.68   9.269

plot(allEffects(CN.RK))

Figure 2
data means and SE for plots
- bar plots of d13C, d15N, TC and TN for all 4 samples types with asterisks for significant effects

## make some means and SE
d13C.m<-aggregate(d13C~Plot+Sample, data=Hak.data, mean)
d15N.m<-aggregate(d15N~Plot+Sample, data=Hak.data, mean)
TN.m<-aggregate(Total.N..mmol.gdw~Plot+Sample, data=Hak.data, mean)
TC.m<-aggregate(Total.C..mmol.gdw~Plot+Sample, data=Hak.data, mean)
CN.m<-aggregate(C.N~Plot+Sample, data=Hak.data, mean)

d13C.SE<-aggregate(d13C~Plot+Sample, data=Hak.data, std.error)
d15N.SE<-aggregate(d15N~Plot+Sample, data=Hak.data, std.error)
TN.SE<-aggregate(Total.N..mmol.gdw~Plot+Sample, data=Hak.data, std.error)
TC.SE<-aggregate(Total.C..mmol.gdw~Plot+Sample, data=Hak.data, std.error)
CN.SE<-aggregate(C.N~Plot+Sample, data=Hak.data, std.error)

d13C.n<-aggregate(d13C~Plot+Sample, data=Hak.data, length)
d15N.n<-aggregate(d15N~Plot+Sample, data=Hak.data, length)

mean.data<- cbind(d13C.m, d15N.m[3], TN.m[3], TC.m[3], d13C.SE[3], d15N.SE[3], TN.SE[3], TC.SE[3], CN.m[3], CN.SE[3])
colnames(mean.data)<-c("Plot", "Sample", "d13C.mean", "d15N.mean", "TN.mean", "TC.mean", "d15N.SE", "d13C.SE", "TN.SE", "TC.SE", "CN.mean", "CN.SE")

# rename Rubus species for plotting
levels(mean.data$Sample)[levels(mean.data$Sample)=="RUBARG"] <- "Rubus spp"
levels(mean.data$Sample)[levels(mean.data$Sample)=="RUBHAW"] <- "Rubus spp"

mean.data$Sample<-factor(mean.data$Sample, levels=c("Soil", "Koa", "Rubus spp"))

### make mean plots
pd <- position_dodge(0.7) #offset for error bars

### d13C
d13C.plot<-ggplot(data=mean.data, aes(x=Sample, y=d13C.mean, fill=Plot, group=Plot))+
  geom_bar(stat="identity", position = pd, width=0.7)+
  geom_errorbar(aes(ymin=d13C.mean-d13C.SE, ymax=d13C.mean+d13C.SE), size=.3, width=0, position=pd) +
  scale_fill_manual(values=c('gray65', "gray89")) +
  xlab("Sample Type") + 
  annotate(geom="text", x=1, y=-27.5, label=expression(paste(bold("*"))), size=4) +
  annotate(geom="text", x=2, y=-33, label=expression(paste(bold("*"))), size=4) +
  annotate(geom="text", x=3, y=-33, label=expression(paste(bold("*"))), size=4) +
  theme(text = element_text(size=8)) +
  ylab(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))) + BW.back2

### d15N
d15N.plot<-ggplot(data=mean.data, aes(x=Sample, y=d15N.mean, fill=Plot, group=Plot))+
  geom_bar(stat="identity", position = pd, width=0.7)+ ylim(c(0,8)) +
  geom_errorbar(aes(ymin=d15N.mean-d15N.SE, ymax=d15N.mean+d15N.SE), size=.3, width=0, position=pd) +
  scale_fill_manual(values=c('gray65', "gray89")) +
  xlab("Sample Type") + 
  theme(text = element_text(size=8)) +
  annotate(geom="text", x=2, y=2.8, label=expression(paste(bold("*"))), size=4) +
  annotate(geom="text", x=3, y=3.3, label=expression(paste(bold("*"))), size=4) +
  ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) + BW.back2

### total Carbon
TC.plot<-ggplot(data=mean.data, aes(x=Sample, y=TC.mean, fill=Plot, group=Plot))+
  geom_bar(stat="identity", position = pd, width=0.7)+ ylim(c(0, 60)) +
  geom_errorbar(aes(ymin=TC.mean-TC.SE, ymax=TC.mean+TC.SE), size=.3, width=0, position=pd) +
  scale_fill_manual(values=c('gray65', "gray89")) +
  xlab("Sample Type") + 
  theme(text = element_text(size=8)) +
  annotate(geom="text", x=1, y=25, label=expression(paste(bold("*"))), size=4) +
  annotate(geom="text", x=3, y=42, label=expression(paste(bold("*"))), size=4) +
  ylab("Total Carbon (mmol/gdw)") + BW.back2

### total Nitrogen
TN.plot<-ggplot(data=mean.data, aes(x=Sample, y=TN.mean, fill=Plot, group=Plot))+
  geom_bar(stat="identity", position = pd, width=0.7)+
  geom_errorbar(aes(ymin=TN.mean-TN.SE, ymax=TN.mean+TN.SE), size=.3, width=0, position=pd) +
  scale_fill_manual(values=c('gray65', "gray89")) +
  xlab("Sample Type") + 
  theme(text = element_text(size=8)) +
  ylab("Total Nitrogen (mmol/gdw)") + BW.back2

# get the legend, # create some space to the left of the legend
bar.legend <- get_legend(
  TN.plot + theme(legend.box.margin = margin(0, 0, 0, 12)) + 
    theme(legend.key.size = unit(0.3, "cm")))

#  figures together
bar.plots<-plot_grid(d15N.plot + theme(legend.position = "none"), 
                     d13C.plot + theme(legend.position = "none"),
                     TN.plot + theme(legend.position = "none"),
                     TC.plot + theme(legend.position = "none"), 
                     nrow=1, ncol=4, labels=c('a', 'b', 'c', 'd'), 
          label_size=8, hjust=-1, vjust= 3)

plot_grid(bar.plots, bar.legend, rel_widths = c(8, 1)) # legend  1/8 size as first obj.

dev.copy(pdf, "figures/Fig 2.mean.data.sample.pdf", width=7, height=3, encod="MacRoman")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

###ISO-KRIG Give the data a “krig-over”.

First, make krigs of the full plot and all data layers

# load data
krig<-read.csv("data/krig.matrix.csv") # matrix of points for isoscape grid

# merge kriging data with isotope data
data.merge<-as.data.frame(join_all(list(Hak.data, krig), by = "Position.point", type='full'))

data.merge<-data.merge[c(-167:-194),] # drop NAs where no samples taken
scape.data<-data.merge

write.csv(scape.data, "output/scape.data.csv")

#################

####d15N isoscape
setup

################# d15N krigs
################# make site maps based on bubble plots

# just AK site
krig.AK<- scape.data %>% filter(Plot=="AK")
krig.AK$Sample<-droplevels(krig.AK$Sample)

# modify one point slightly due to overlap of coordinate system
krig.AK$y.matrix[24]=2.4; krig.AK$y.meter[24]=2.4
krig.AK$x.matrix[24]=5.4; krig.AK$x.meter[24]=5.4

AK.map<-krig.AK %>% as.data.frame %>% 
  ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d15N, color=Sample), alpha=3/4) + 
  ggtitle("AK plot--d15N (permil)") + coord_equal() + theme_bw()

# just RK site
krig.RK<- scape.data %>% filter(Plot=="RK")
krig.RK$Sample<-droplevels(krig.RK$Sample)

RK.map<-krig.RK %>% as.data.frame %>% 
  ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d15N, color=Sample), alpha=3/4) +
  ggtitle("RK plot--d15N (permil)") + coord_equal() + theme_bw()

plot_grid(RK.map, AK.map, ncol=2, nrow=1)

dev.copy(pdf, "figures/dot.scape.pdf", width=9, height=4, encod="MacRoman")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

#####################
#####################

AK d15N isoscape

## Krig AK sites
coordinates(krig.AK)<- ~x.meter + y.meter

# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points

## expand binding box 
#krig.AK@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)

##   min max
## x  -1  36
## y  -1  21

krig.AK@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner

# autokrige
d15N.AK.auto <- autoKrige(d15N ~ 1, krig.AK, Grid.AK) # ordinary kriging

## [using ordinary kriging]

plot(d15N.AK.auto, sp.layout = list(pts = list("sp.points", krig.AK)))

#d15N.AK.auto$krige_output

# make to dataframes for lm, combine
AK.pred.N <- d15N.AK.auto[1] %>% as.data.frame() # model predictions
AK.df.N <- krig.AK %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.AK.df.N<- left_join(AK.df.N, AK.pred.N, by = c("x.meter"="krige_output.x")) # longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.df.N)
summary(mod) #R squared = 0.21

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.AK.df.N)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6314 -0.4768  0.0794  0.4671  4.0753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.818896   0.009430  404.99   <2e-16 ***
## d15N        0.184962   0.001937   95.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8717 on 35278 degrees of freedom
## Multiple R-squared:  0.2054, Adjusted R-squared:  0.2054 
## F-statistic:  9119 on 1 and 35278 DF,  p-value: < 2.2e-16

mean(pred.AK.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 4.4 d15N

## [1] 4.602741

mean(pred.AK.df.N$krige_output.var1.stdev, na.rm=T) # 2.1

## [1] 2.088436

# inspect model
# plot(pred.AK.df.N$krige_output.var1.pred ~ pred.AK.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

### map in automap
#automapPlot(d15N.AK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.AK, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# note that the prediction "var1.pred" component of the krig is a SpatialPixelsDataFrame. Need to call the krig, then the component df("krige_output"), then the output column ("var1.pred") to plot in spplot

# variance
# spplot(d15N.AK.auto$krige_output,"var1.stdev")

#predicted krige
my.palette2 <- brewer.pal(n = 6, name = "BuGn")
my.palette3 <- colorRampPalette(c("firebrick1", "darkseagreen1", "azure1"))(4)
col.scheme.N <- colorRampPalette(my.palette3)


plotd15N.krig.AK<-spplot(d15N.AK.auto$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", krig.AK, pch=c(16,17,1)[krig.AK$Sample], 
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)


###########################

RK d15N isoscape

###########################
# RK Site Krig
###########################
## Krig RK site 
krig.RK<- scape.data %>% filter(Plot=="RK")
krig.RK$Sample<-droplevels(krig.RK$Sample)


coordinates(krig.RK)<- ~x.meter + y.meter

Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points

## expand binding box 
#krig.RK@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)

##   min max
## x  -1  21
## y  -1  36

Grid.RK@bbox<-bbox(S) # expand for plot corner

# autokrige
# 'fix.values' = The items describe the fixed value for the nugget (y-intercept), range (of interprolation) and sill(asymptote) respectively. Setting the value to NA means that the value is not fixed. Is passed on to autofitVariogram.
# The nugget is the y-intercept of the variogram. In practical terms, the nugget represents the small-scale variability of the data. A portion of that short range variability can be the result of measurement error. 
# The range is the distance after which the variogram levels off. The physical meaning of the range is that pairs of points that are this distance or greater apart are not spatially correlated. 
# The sill is the total variance contribution, or the maximum variability between pairs of points. 

# increasing the range gave better fit than auto
d15N.RK.auto <- autoKrige(d15N ~ 1, krig.RK, Grid.RK, fix.value=c(NA, 2, NA)) # ordinary kriging

## [using ordinary kriging]

plot(d15N.RK.auto, sp.layout = list(pts = list("sp.points", krig.RK)))

#d15N.RK.auto$krige_output

# make to dataframes for lm, combine
RK.pred.N <- d15N.RK.auto[1] %>% as.data.frame() # model predictions
RK.df.N <- krig.RK %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
# longest axis here is y
pred.RK.df.N<- left_join(RK.df.N, RK.pred.N, by = c("x.meter"="krige_output.x"))
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.df.N)
summary(mod) #R squared = 0.21

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.RK.df.N)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7946 -0.4047  0.0087  0.3657  3.4740 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.026386   0.007487   537.8   <2e-16 ***
## d15N        0.184935   0.001434   128.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7525 on 62983 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2089, Adjusted R-squared:  0.2088 
## F-statistic: 1.663e+04 on 1 and 62983 DF,  p-value: < 2.2e-16

mean(pred.RK.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 5.0 d15N

## [1] 4.911052

mean(pred.RK.df.N$krige_output.var1.stdev, na.rm=T) # 1.9

## [1] 1.90181

# inspect model
# plot(pred.RK.df.N$krige_output.var1.pred ~ pred.RK.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

### map in automap
#automapPlot(d15N.RK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.RK), col.regions=my.palette)

# variance
# spplot(d15N.RK.auto$krige_output,"var1.stdev")

#predicted krige
plotd15N.krig.RK<-spplot(d15N.RK.auto$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", krig.RK, pch=c(16,17,1)[krig.RK$Sample], 
                        col="black", cex=0.5, alpha=0.5), 
       main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.7))


#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
  print(x)
}
### combined Ak-Rk krig d15N plot 
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK, plotd15N.krig.RK, ncol=2, nrow=1, widths = c(1.2,1))

dev.copy(png, "figures/RKAK.d15N.krig3.png", width = 8, height = 4, units = 'in', res = 300)

## quartz_off_screen 
##                 3

dev.off()

## quartz_off_screen 
##                 2

d13C isoscape

setup

################# d13C krigs
################# make site maps based on bubble plots

AK.map<-krig.AK %>% as.data.frame %>% 
  ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d13C, color=Sample), alpha=3/4) + 
  ggtitle("AK plot--d13C (permil)") + coord_equal() + theme_bw()


RK.map<-krig.RK %>% as.data.frame %>% 
  ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d13C, color=Sample), alpha=3/4) +
  ggtitle("RK plot--d13C (permil)") + coord_equal() + theme_bw()

plot_grid(RK.map, AK.map, ncol=2, nrow=1)

dev.copy(pdf, "figures/dot.scape.d13C.pdf", width=9, height=4, encod="MacRoman")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

#####################
#####################

AK d13C isoscape

# autokrige
d13C.AK.auto <- autoKrige(d13C ~ 1, krig.AK, Grid.AK) # ordinary kriging

## [using ordinary kriging]

plot(d13C.AK.auto, sp.layout = list(pts = list("sp.points", krig.AK)))

#d13C.AK.auto$krige_output

# make to dataframes for lm, combine
AK.pred.C <- d13C.AK.auto[1] %>% as.data.frame() # model predictions
AK.df.C <- krig.AK %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.AK.df.C<- left_join(AK.df.C, AK.pred.C, by = c("x.meter"="krige_output.x")) # longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.AK.df.C)
summary(mod) #R squared = 0.24

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.AK.df.C)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3830 -0.4285  0.0828  0.5165  4.2039 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -21.5295     0.0423  -509.0   <2e-16 ***
## d13C          0.1677     0.0016   104.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7706 on 35278 degrees of freedom
## Multiple R-squared:  0.2375, Adjusted R-squared:  0.2374 
## F-statistic: 1.099e+04 on 1 and 35278 DF,  p-value: < 2.2e-16

mean(pred.AK.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -25.9 d13C

## [1] -25.9423

mean(pred.AK.df.C$krige_output.var1.stdev, na.rm=T) # 2.3

## [1] 2.34059

# inspect model
# plot(pred.AK.df.C$krige_output.var1.pred ~ pred.AK.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d13C.AK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.AK, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d13C.AK.auto$krige_output,"var1.stdev")

#predicted krige
my.palette4 <- colorRampPalette(c("firebrick1", "paleturquoise2", "azure1"))(4)
col.scheme.C <- colorRampPalette(my.palette4)


plotd13C.krig.AK<-spplot(d13C.AK.auto$krige_output["var1.pred"], col.regions=col.scheme.C, 
       sp.layout = list("sp.points", krig.AK, pch=c(16,17,1)[krig.AK$Sample],
                        col="black", cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{13}, C, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)


###########################

RK d13C isoscape

###########################
# RK Site Krig
###########################
## Krig RK site 

# autokrige
d13C.RK.auto <- autoKrige(d13C ~ 1, krig.RK, Grid.RK, fix.value=c(NA, 1, NA)) # ordinary kriging

## [using ordinary kriging]

plot(d13C.RK.auto, sp.layout = list(pts = list("sp.points", krig.RK)))

#d13C.RK.auto$krige_output

# make to dataframes for lm, combine
RK.pred.C <- d13C.RK.auto[1] %>% as.data.frame() # model predictions
RK.df.C <- krig.RK %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
# longest axis here is y
pred.RK.df.C<- left_join(RK.df.C, RK.pred.C, by = c("x.meter"="krige_output.x"))
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.RK.df.C)
summary(mod) #R squared = 0.33

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.RK.df.C)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7059 -0.7070  0.0331  0.7676  4.2973 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -20.013291   0.043240  -462.8   <2e-16 ***
## d13C          0.271933   0.001551   175.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.093 on 62983 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3279, Adjusted R-squared:  0.3279 
## F-statistic: 3.073e+04 on 1 and 62983 DF,  p-value: < 2.2e-16

mean(pred.RK.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -27.6 d13C

## [1] -27.55484

mean(pred.RK.df.C$krige_output.var1.stdev, na.rm=T) # 2.8

## [1] 2.50984

# inspect model
# plot(pred.RK.df.C$krige_output.var1.pred ~ pred.RK.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# map in automap
#automapPlot(d13C.RK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.RK), col.regions=my.palette)

# variance
# spplot(d13C.RK.auto$krige_output,"var1.stdev")

#predicted krige
plotd13C.krig.RK<-spplot(d13C.RK.auto$krige_output["var1.pred"], col.regions=col.scheme.C, 
       sp.layout = list("sp.points", krig.RK, pch=c(16,17,1)[krig.RK$Sample],
                        col="black", cex=0.5, alpha=0.5), 
       main=list(label=expression(paste(delta^{13}, C, " - RK (Remnant Forest)")), cex=0.7))



#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
  print(x)
}
### combined Ak-Rk krig d13C plot 
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd13C.krig.AK,plotd13C.krig.RK, ncol=2, nrow=1, widths = c(1.2,1))

dev.copy(png, "figures/RKAK.d13C.krig2.png", width = 8, height = 4, units = 'in', res = 300)

## quartz_off_screen 
##                 3

dev.off()

## quartz_off_screen 
##                 2

Density plots: all samples

Density plots using predictions of N and C values for all sample types

###################
# make density plots using predictions
# combine the predictions for each plot
AK.df.out.N<-as.data.frame(AK.pred.N$krige_output.var1.pred); AK.df.out.N$Plot<-"AK"
RK.df.out.N<-as.data.frame(RK.pred.N$krige_output.var1.pred); RK.df.out.N$Plot<-"RK"; 
colnames(AK.df.out.N)<-c("d15N.pred", "Plot")
colnames(RK.df.out.N)<-c("d15N.pred", "Plot")

# new dataframe
pred.dat.N<-rbind(AK.df.out.N, RK.df.out.N)
pred.dat.N$Plot<-as.factor(pred.dat.N$Plot)
pred.dat.N %>% 
  group_by(Plot) %>%
  summarise(no_rows = length(Plot))

##   no_rows
## 1  653562

#man whitney for significance (not normal data)
mwu(pred.dat.N, d15N.pred, Plot) # signif difference in predictions (p<0.001)

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 326781) | 2 = RK (n = 326781):
##   U = 83294984480.000, W = 29901910109.000, p < 0.001, Z = -308.086
##   effect-size r =      0.381
##    rank-mean(1) = 254895.43
##    rank-mean(2) = 398667.57

# Density plot

#hist(pred.dat.N$d15N.pred[(pred.dat.N$Plot=="AK")], ylim=c(0,150000), xlab="d15N-predict")
#hist(pred.dat.N$d15N.pred[(pred.dat.N$Plot=="RK")], ylim=c(0,150000), xlab="d15N-predict")

# plot.means
predict.mean.N <- ddply(pred.dat.N, "Plot", summarise, grp.mean=mean(d15N.pred, na.rm=TRUE))

density.predict.N<-ggplot(pred.dat.N, aes(x=d15N.pred)) +
  geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+
  xlab(expression(paste(delta^{15}, N["predicted"], sp="")))+
  scale_color_manual(values=c("#E69F00", "#56B4E9")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
  geom_vline(data=predict.mean.N, aes(xintercept=grp.mean, color=Plot),
             linetype="dashed", lwd=0.2) + BW.back2


######## Carbon


###################
# make density plots using predictions
# combine the predictions for each plot
AK.df.out.C<-as.data.frame(AK.pred.C$krige_output.var1.pred); AK.df.out.C$Plot<-"AK"
RK.df.out.C<-as.data.frame(RK.pred.C$krige_output.var1.pred); RK.df.out.C$Plot<-"RK"; 
colnames(AK.df.out.C)<-c("d13C.pred", "Plot")
colnames(RK.df.out.C)<-c("d13C.pred", "Plot")

# new dataframe
pred.dat.C<-rbind(AK.df.out.C, RK.df.out.C)
pred.dat.C$Plot<-as.factor(pred.dat.C$Plot)
pred.dat.C %>% 
  group_by(Plot) %>%
  summarise(no_rows = length(Plot))

##   no_rows
## 1  653562

#man whitney for significance (not normal data)
mwu(pred.dat.C, d13C.pred, Plot) # signif difference in predictions (p<0.001)

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 326781) | 2 = RK (n = 326781):
##   U = 147098098591.500, W = 93705024220.500, p < 0.001, Z = 528.598
##   effect-size r =      0.654
##    rank-mean(1) = 450142.75
##    rank-mean(2) = 203420.25

# Density plot

#hist(pred.dat.C$d13C.pred[(pred.dat.C$Plot=="AK")], ylim=c(0,150000), xlab="d13C-predict")
#hist(pred.dat.C$d13C.pred[(pred.dat.C$Plot=="RK")], ylim=c(0,150000), xlab="d13C-predict")

# plot.means
predict.mean.C <- ddply(pred.dat.C, "Plot", summarise, grp.mean=mean(d13C.pred, na.rm=TRUE))

density.predict.C<-ggplot(pred.dat.C, aes(x=d13C.pred)) +
  geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+
  xlab(expression(paste(delta^{13}, C["predicted"], sp="")))+
  scale_color_manual(values=c("#E69F00", "#56B4E9")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
  geom_vline(data=predict.mean.C, aes(xintercept=grp.mean, color=Plot),
             linetype="dashed", lwd=0.2) + BW.back2

# get the legend, # create some space to the left of the legend
density.legend <- get_legend(
  density.predict.C + theme(legend.box.margin = margin(0, 0, 0, 12)) + 
    theme(legend.key.size = unit(0.3, "cm")))

gridExtra::grid.arrange(density.predict.N + theme(legend.position = "none"),
                        density.predict.C + theme(legend.position = "none"),
                        density.legend,
                        ncol=3, nrow=1, widths = c(1,1,0.4))

dev.copy(pdf, "figures/Fig S2ab.mod.pred.NC.pdf", width = 6, height = 2.5)

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

Soil krigs

d15N Soil isoscape

Subset data and only examine SOIL.
- Restored forest (AK) and Remnant Forest (RK)

###### AK site d15N

## d15N SOIL plot
krig.AK<- scape.data %>% filter(Plot=="AK") # just AK
AK.soil<-  krig.AK[(krig.AK$Sample=="Soil"),] #just soil
AK.soil$Sample<-droplevels(AK.soil$Sample)

## Krig AK sites
coordinates(AK.soil)<- ~x.meter + y.meter

# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points

## expand binding box 
#AK.soil@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)

##   min max
## x  -1  36
## y  -1  21

AK.soil@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner

# autokrige
d15N.AK.soil <- autoKrige(d15N ~ 1, AK.soil, Grid.AK) # ordinary kriging

## [using ordinary kriging]

plot(d15N.AK.soil, sp.layout = list(pts = list("sp.points", AK.soil)))

#d15N.AK.soil$krige_output

# make to dataframes for lm, combine
AK.soil.pred.N <- d15N.AK.soil[1] %>% as.data.frame() # model predictions
AK.soil.df.N <- AK.soil %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.AK.soil.df.N<- left_join(AK.soil.df.N, AK.soil.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.soil.df.N)
summary(mod) #R squared = 0.015

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.AK.soil.df.N)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3313 -0.5382 -0.0558  0.6719  3.0052 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.31500    0.03730  142.49   <2e-16 ***
## d15N         0.11105    0.00609   18.24   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9752 on 21166 degrees of freedom
## Multiple R-squared:  0.01547,    Adjusted R-squared:  0.01542 
## F-statistic: 332.6 on 1 and 21166 DF,  p-value: < 2.2e-16

mean(pred.AK.soil.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 6.0 d15N

## [1] 5.984172

mean(pred.AK.soil.df.N$krige_output.var1.stdev, na.rm=T) # 0.35

## [1] 0.3514554

# inspect model
# plot(pred.AK.soil.df.N$krige_output.var1.pred ~ pred.AK.soil.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d15N.AK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", AK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d15N.AK.soil$krige_output,"var1.stdev")

#predicted krige
plotd15N.krig.AK.soil<-spplot(d15N.AK.soil$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", AK.soil, pch=1, 
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)


###########################
# RK Site Krig
###########################

## d15N SOIL plot
krig.RK<- scape.data %>% filter(Plot=="RK") # just AK
RK.soil<-  krig.RK[(krig.RK$Sample=="Soil"),] #just soil
RK.soil$Sample<-droplevels(RK.soil$Sample)

## Krig RK sites
coordinates(RK.soil)<- ~x.meter + y.meter

# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points

## expand binding box 
#RK.soil@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)

##   min max
## x  -1  21
## y  -1  36

RK.soil@bbox<-bbox(S) # expanded binding box for data
Grid.RK@bbox<-bbox(S) # expand for plot corner

# autokrige
d15N.RK.soil <- autoKrige(d15N ~ 1, RK.soil, Grid.RK) # ordinary kriging

## [using ordinary kriging]

plot(d15N.RK.soil, sp.layout = list(pts = list("sp.points", RK.soil)))

#d15N.RK.soil$krige_output

# make to dataframes for lm, combine
RK.soil.pred.N <- d15N.RK.soil[1] %>% as.data.frame() # model predictions
RK.soil.df.N <- RK.soil %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.RK.soil.df.N<- left_join(RK.soil.df.N, RK.soil.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.soil.df.N)
summary(mod) #R squared = 0.01

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.RK.soil.df.N)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.87178 -0.34756  0.02653  0.50770  1.97596 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.124328   0.022592  271.08   <2e-16 ***
## d15N        0.052800   0.003492   15.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.708 on 35566 degrees of freedom
## Multiple R-squared:  0.006386,   Adjusted R-squared:  0.006358 
## F-statistic: 228.6 on 1 and 35566 DF,  p-value: < 2.2e-16

mean(pred.RK.soil.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 6.4 d15N

## [1] 6.461148

mean(pred.RK.soil.df.N$krige_output.var1.stdev, na.rm=T) # 0.8

## [1] 0.7599088

# inspect model
# plot(pred.RK.soil.df.N$krige_output.var1.pred ~ pred.RK.soil.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d15N.RK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", RK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d15N.RK.soil$krige_output,"var1.stdev")

#predicted krige
plotd15N.krig.RK.soil<-spplot(d15N.RK.soil$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", RK.soil, pch=1, 
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)


#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
  print(x)
}
### combined Ak-Rk krig d15N plot 
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK.soil, plotd15N.krig.RK.soil, ncol=2, nrow=1, widths = c(1.2,1))

dev.copy(png, "figures/RKAK.d15N.soils.png", width = 8, height = 4, units = 'in', res = 300)

## quartz_off_screen 
##                 3

dev.off()

## quartz_off_screen 
##                 2

d13C Soil isoscape

###### AK site d13C

# autokrige
d13C.AK.soil <- autoKrige(d13C ~ 1, AK.soil, Grid.AK) # ordinary kriging

## [using ordinary kriging]

plot(d13C.AK.soil, sp.layout = list(pts = list("sp.points", AK.soil)))

#d13C.AK.soil$krige_output

# make to dataframes for lm, combine
AK.soil.pred.C <- d13C.AK.soil[1] %>% as.data.frame() # model predictions
AK.soil.df.C <- AK.soil %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.AK.soil.df.C<- left_join(AK.soil.df.C, AK.soil.pred.C, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.AK.soil.df.C)
summary(mod) #R squared = 0.008

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.AK.soil.df.C)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.61780 -0.52385  0.02572  0.48253  2.96360 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -22.734708   0.130976 -173.58   <2e-16 ***
## d13C          0.069928   0.005368   13.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8341 on 21166 degrees of freedom
## Multiple R-squared:  0.007955,   Adjusted R-squared:  0.007908 
## F-statistic: 169.7 on 1 and 21166 DF,  p-value: < 2.2e-16

mean(pred.AK.soil.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -24.4 d13C

## [1] -24.43935

mean(pred.AK.soil.df.C$krige_output.var1.stdev, na.rm=T) # 0.6

## [1] 0.6411003

# inspect model
# plot(pred.AK.soil.df.C$krige_output.var1.pred ~ pred.AK.soil.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d13C.AK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", AK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d13C.AK.soil$krige_output,"var1.stdev")

#predicted krige
plotd13C.krig.AK.soil<-spplot(d13C.AK.soil$krige_output["var1.pred"], col.regions=col.scheme.C, 
       sp.layout = list("sp.points", AK.soil, pch=1, 
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{13}, C, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)



###########################
# RK Site Krig
###########################

# autokrige
d13C.RK.soil <- autoKrige(d13C ~ 1, RK.soil, Grid.RK) # ordinary kriging

## [using ordinary kriging]

plot(d13C.RK.soil, sp.layout = list(pts = list("sp.points", RK.soil)))

#d13C.RK.soil$krige_output

# make to dataframes for lm, combine
RK.soil.pred.C <- d13C.RK.soil[1] %>% as.data.frame() # model predictions
RK.soil.df.C <- RK.soil %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.RK.soil.df.C<- left_join(RK.soil.df.C, RK.soil.pred.C, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.RK.soil.df.C)
summary(mod) #R squared = 0.16

## 
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.RK.soil.df.C)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.87389 -0.51641 -0.01577  0.55010  2.42974 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -22.840353   0.093601 -244.02   <2e-16 ***
## d13C          0.097864   0.003685   26.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6794 on 35566 degrees of freedom
## Multiple R-squared:  0.01944,    Adjusted R-squared:  0.01941 
## F-statistic: 705.2 on 1 and 35566 DF,  p-value: < 2.2e-16

mean(pred.RK.soil.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions 4.4 d13C

## [1] -25.32413

mean(pred.RK.soil.df.C$krige_output.var1.stdev, na.rm=T) # 2.1

## [1] 0.6704758

# inspect model
# plot(pred.RK.soil.df.C$krige_output.var1.pred ~ pred.RK.soil.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d13C.RK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", RK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d13C.RK.soil$krige_output,"var1.stdev")

#predicted krige
plotd13C.krig.RK.soil<-spplot(d13C.RK.soil$krige_output["var1.pred"], col.regions=col.scheme.C, 
       sp.layout = list("sp.points", RK.soil, pch=1, 
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{13}, C, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)


#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
  print(x)
}
### combined Ak-Rk krig d13C plot 
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd13C.krig.AK.soil, plotd13C.krig.RK.soil, ncol=2, nrow=1, widths = c(1.2,1))

dev.copy(png, "figures/RKAK.d13C.soils.png", width = 8, height = 4, units = 'in', res = 300)

## quartz_off_screen 
##                 3

dev.off()

## quartz_off_screen 
##                 2

Density plots for soil isoscapes

# make density plots using predictions
# combine the predictions for each plot


############## Nitrogen plot
AK.df.soil.N<-as.data.frame(pred.AK.soil.df.N$krige_output.var1.pred); AK.df.soil.N$Plot<-"AK"
RK.df.soil.N<-as.data.frame(pred.RK.soil.df.N$krige_output.var1.pred); RK.df.soil.N$Plot<-"RK" 
colnames(AK.df.soil.N)<-c("d15N.pred", "Plot")
colnames(RK.df.soil.N)<-c("d15N.pred", "Plot")

# new dataframe
pred.soil.N<-rbind(AK.df.soil.N, RK.df.soil.N)
pred.soil.N$Plot<-as.factor(pred.soil.N$Plot)
pred.soil.N %>% 
  group_by(Plot) %>%
  summarise(no_rows = length(Plot))

##   no_rows
## 1   56736

#man whitney for significance (not normal data)
mwu(pred.soil.N, d15N.pred, Plot) # signif difference in predictions (p<0.001)

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 21168) | 2 = RK (n = 35568):
##   U = 474341208.000, W = 250288512.000, p < 0.001, Z = -66.868
##   effect-size r =     0.281
##    rank-mean(1) = 22408.41
##    rank-mean(2) = 31915.60

# Density plot

#hist(pred.soil.N$d15N.pred[(pred.soil.N$Plot=="AK")], ylim=c(0,10000), xlab="d15N-predict")
#hist(pred.soil.N$d15N.pred[(pred.soil.N$Plot=="RK")], ylim=c(0,10000), xlab="d15N-predict")

# plot.means
predict.mean.N <- ddply(pred.soil.N, "Plot", summarise, grp.mean=mean(d15N.pred, na.rm=TRUE))

density.predict.Nsoil<-ggplot(pred.soil.N, aes(x=d15N.pred)) +
  geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+ ylim(0,1.2) +
  xlab(expression(paste(delta^{15}, N["predicted"], sp="")))+
  scale_color_manual(values=c("#E69F00", "#56B4E9")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
  geom_vline(data=predict.mean.N, aes(xintercept=grp.mean, color=Plot),
             linetype="dashed", lwd=0.2) + BW.back2


############## Carbon
AK.df.soil.C<-as.data.frame(pred.AK.soil.df.C$krige_output.var1.pred); AK.df.soil.C$Plot<-"AK"
RK.df.soil.C<-as.data.frame(pred.RK.soil.df.C$krige_output.var1.pred); RK.df.soil.C$Plot<-"RK" 
colnames(AK.df.soil.C)<-c("d13C.pred", "Plot")
colnames(RK.df.soil.C)<-c("d13C.pred", "Plot")

# new dataframe
pred.soil.C<-rbind(AK.df.soil.C, RK.df.soil.C)
pred.soil.C$Plot<-as.factor(pred.soil.C$Plot)
pred.soil.C %>% 
  group_by(Plot) %>%
  summarise(no_rows = length(Plot))

##   no_rows
## 1   56736

#man whitney for significance (not normal data)
mwu(pred.soil.C, d13C.pred, Plot) # signif difference in predictions (p<0.001)

## 
## # Mann-Whitney-U-Test
## 
## Groups 1 = AK (n = 21168) | 2 = RK (n = 35568):
##   U = 819376896.000, W = 595324200.000, p < 0.001, Z = 116.006
##   effect-size r =     0.487
##    rank-mean(1) = 38708.28
##    rank-mean(2) = 22214.87

# Density plot

#hist(pred.soil.C$d13C.pred[(pred.soil.C$Plot=="AK")], ylim=c(0,10000), xlab="d13C-predict")
#hist(pred.soil.C$d13C.pred[(pred.soil.C$Plot=="RK")], ylim=c(0,10000), xlab="d13C-predict")

# plot.means
predict.mean.C <- ddply(pred.soil.C, "Plot", summarise, grp.mean=mean(d13C.pred, na.rm=TRUE))

density.predict.Csoil<-ggplot(pred.soil.C, aes(x=d13C.pred)) +
  geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+ ylim(0,0.8) +
  xlab(expression(paste(delta^{13}, C["predicted"], sp="")))+
  scale_color_manual(values=c("#E69F00", "#56B4E9")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
  geom_vline(data=predict.mean.C, aes(xintercept=grp.mean, color=Plot),
             linetype="dashed", lwd=0.2) + BW.back2

# get the legend, # create some space to the left of the legend
density.legend <- get_legend(
  density.predict.Nsoil + theme(legend.box.margin = margin(0, 0, 0, 12)) + 
    theme(legend.key.size = unit(0.3, "cm")))

gridExtra::grid.arrange(density.predict.Nsoil + theme(legend.position = "none"),
                        density.predict.Csoil + theme(legend.position = "none"),
                        density.legend,
                        ncol=3, nrow=1, widths = c(1,1,0.3))
dev.copy(pdf, "figures/Fig Sx.mod.pred.pdf", width = 5, height = 2.5)

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

#### d15N Plants isoscape
Subset data and only examine PLANTS     
- *Restored forest (AK)* and *Remnant Forest (RK)*
###### AK site d15N

## d15N plants plot
krig.AK<- scape.data %>% filter(Plot=="AK") # just AK
AK.plants<-  krig.AK[!(krig.AK$Sample=="Soil"),] #just plants
AK.plants$Sample<-droplevels(AK.plants$Sample)

## Krig AK sites
coordinates(AK.plants)<- ~x.meter + y.meter

# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points

## expand binding box 
#AK.plants@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
AK.plants@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner

# autokrige
d15N.AK.plants <- autoKrige(d15N ~ 1, AK.plants, Grid.AK, fix.value=c(0.2, 3, 0.6)) # ordinary kriging
plot(d15N.AK.plants, sp.layout = list(pts = list("sp.points", AK.plants)))
#d15N.AK.plants$krige_output

# make to dataframes for lm, combine
AK.plants.pred.N <- d15N.AK.plants[1] %>% as.data.frame() # model predictions
AK.plants.df.N <- AK.plants %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.AK.plants.df.N<- left_join(AK.plants.df.N, AK.plants.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.plants.df.N)
summary(mod) #R squared = 0.19

mean(pred.AK.plants.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 1.6 d15N
mean(pred.AK.plants.df.N$krige_output.var1.stdev, na.rm=T) # 0.7

# inspect model
# plot(pred.AK.plants.df.N$krige_output.var1.pred ~ pred.AK.plants.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d15N.AK.plants$krige_output, "var1.pred", sp.layout = list("sp.points", AK.plants, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d15N.AK.plants$krige_output,"var1.stdev")

#predicted krige
plotd15N.krig.AK.plants<-spplot(d15N.AK.plants$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", AK.plants, pch=c(16,17)[AK.plants$Sample],
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)

plotd15N.krig.AK.plants


###########################
# RK Site Krig
###########################

## d15N plants plot
krig.RK<- scape.data %>% filter(Plot=="RK") # just AK
RK.plants<-  krig.RK[!(krig.RK$Sample=="Soil"),] #just plants
RK.plants$Sample<-droplevels(RK.plants$Sample)

## Krig RK sites
coordinates(RK.plants)<- ~x.meter + y.meter

# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points

## expand binding box 
#RK.plants@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
RK.plants@bbox<-bbox(S) # expanded binding box for data
Grid.RK@bbox<-bbox(S) # expand for plot corner

# autokrige
d15N.RK.plants <- autoKrige(d15N ~ 1, RK.plants, Grid.RK) # ordinary kriging
plot(d15N.RK.plants, sp.layout = list(pts = list("sp.points", RK.plants)))
#d15N.RK.plants$krige_output

# make to dataframes for lm, combine
RK.plants.pred.N <- d15N.RK.plants[1] %>% as.data.frame() # model predictions
RK.plants.df.N <- RK.plants %>% as.data.frame() # dataframe

# make dataframe with data and predictions (krig) to assess fit
pred.RK.plants.df.N<- left_join(RK.plants.df.N, RK.plants.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.plants.df.N)
summary(mod) #R squared = 0.15

mean(pred.RK.plants.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 2.8 d15N
mean(pred.RK.plants.df.N$krige_output.var1.stdev, na.rm=T) # 0.5

# inspect model
# plot(pred.RK.plants.df.N$krige_output.var1.pred ~ pred.RK.plants.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")

# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")

# map in automap
#automapPlot(d15N.RK.plants$krige_output, "var1.pred", sp.layout = list("sp.points", RK.plants, pch=16, col="gray20", cex=0.5), col.regions=my.palette)

# map in ssplot
# variance
# spplot(d15N.RK.plants$krige_output,"var1.stdev")

#predicted krige
plotd15N.krig.RK.plants<-spplot(d15N.RK.plants$krige_output["var1.pred"], col.regions=col.scheme.N, 
       sp.layout = list("sp.points", RK.plants, pch=c(16,17)[RK.plants$Sample],
                        col=c("black"), cex=0.5, alpha=0.5),
       main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)
plotd15N.krig.RK.plants


#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
  print(x)
}
### combined Ak-Rk krig d15N plot 
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK.plants, plotd15N.krig.RK.plants, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d15N.plants.png", width = 8, height = 4, units = 'in', res = 300)
dev.off()

#####################
# Carbon kriging, leave alone for now.
#####################

## Krig AK site d13C
coordinates(krig.AK)<- ~x.meter + y.meter
col.scheme.C <- colorRampPalette(c('coral', 'cadetblue2', 'gray5'))(5)

Grid.AK <- spsample(krig.AK, type='regular', n=1e5)
gridded(Grid.AK) <- TRUE

## d13C plot
d13C.krig.AK <- krige(d13C ~ 1, krig.AK, Grid.AK)
plotd13C.krig.AK<-spplot(d13C.krig.AK["var1.pred"], col.regions=colorRampPalette(col.scheme.C), 
       main=expression(paste(delta^{13}, C, "--AK (Afforested Koa Plot)")))

## Krig RK site d13C
coordinates(krig.RK)<- ~x.meter + y.meter
GridRK <- spsample(krig.RK, type='regular', n=1e5)
gridded(GridRK) <- TRUE

## d13C plot
d13C.krig.RK <- krige(d13C ~ 1, krig.RK, GridRK)
plotd13C.krig.RK<-spplot(d13C.krig.RK["var1.pred"], col.regions=colorRampPalette(col.scheme.C), 
       main=expression(paste(delta^{13}, C, "--RK (Remnant Koa Plot)")))